!=======Programa que faz a trajetória do chute do Roberto Carlos========

PROGRAM exer

!=======================================================================

IMPLICIT NONE

!=======================Declaração de variáveis=========================

      INTEGER :: I
      REAL(8) :: betamb, tetha0, phi0, vi, gammamb 
      REAL(8), DIMENSION(3) :: veli1, posi1, posi, veli
      REAL(8), PARAMETER :: v0=(100.0d0/3.6d0), deltat=0.01d0, g=9.8d0
      REAL(8), PARAMETER :: a1=0.0039d0, a2=0.0058d0, vd=35.0d0
      REAL(8), PARAMETER :: delta=5.0d0, omega=(39.0d0*2.0d0*3.1415926d0)

!============Abertura de arquivos e atribuições de variáveis============

      WRITE(*,*) "Entre com o parametro beta0/mb:"
      READ(*,*) betamb
      WRITE(*,*) "Entre com o ângulo inicial tetha0:"
      READ(*,*) tetha0
      WRITE(*,*) "Entre com o ângulo inicial phi0:"
      READ(*,*) phi0

      OPEN(UNIT=10, FILE="chute_out.dat", STATUS="unknown")

      posi = 0.0d0
      posi1 = 0.0d0
      veli = 0.0d0

!=========================Realiza as iterações==========================

      veli = (/ v0*DSIN(tetha0)*DCOS(phi0), v0*DSIN(tetha0)*DSIN(phi0),&
      &v0*DCOS(tetha0) /)
      DO WHILE (posi1(1) <= 40.0d0)
            !Para o plot do 3d
            !WRITE(10,*) posi1
            !Para o plot do 2d
            WRITE(10,*) posi1(1), posi1(2)
            vi = SQRT( veli(1)**2 + veli(2)**2 + veli(3)**2 )
            gammamb = a1 + (a2/(1.0d0+DEXP((vi-vd)/delta)))
            !Itera velocidade em x
            veli1(1) = veli(1) - (gammamb*vi*veli(1)+betamb*omega*&
            &veli(2))*deltat
            !Itera velocidade em y
            veli1(2) = veli(2) - (gammamb*vi*veli(2)-betamb*omega*&
            veli(1))*deltat
            !Itera velocidade em z
            veli1(3) = veli(3) - (g+gammamb*vi*veli(3))*deltat
            !Calcula a posição x
            posi1(1) = posi(1) + veli(1)*deltat
            !Calcula a posição y
            posi1(2) = posi(2) + veli(2)*deltat
            !Calcula a posição z
            posi1(3) = posi(3) + veli(3)*deltat
            !Testa para ver se z é negativo, se for sai do loop
            IF (posi(3) < 0.0) exit
            !Prepara variáveis para proxima iteração
            posi = posi1
            veli = veli1
      END DO
      !Para o plot do 3d
      !WRITE(10,*) posi
      !Para o plot do 2d
      WRITE(10,*) posi(1), posi(2)

!===================Testa para ver se o gol foi feito===================

      IF (posi(1) > 40.0d0 .AND. posi(3) < 2.5d0 .AND. posi(2) > 4.0d0&
      & .AND. posi(2) < 10.0d0) THEN
            WRITE(*,*) "Você fez o gol? R: sim"
      ELSE
            WRITE(*,*) "Você fez o gol? R: nao"
      END IF

!=======================================================================

END PROGRAM exer

!=======================================================================
